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Abstract Applied problems of oil and gas recovery are studied numerically 
using the mathematical models of multiphase fluid flows in porous media. 
f~^ ■ The basic model includes the continuity equations and the Darcy laws for 

Psl \ each phase, as well as the algebraic expression for the sum of saturations. 

Primary computational algorithms are implemented for such problems using 
the pressure equation. In this paper, we highlight the basic properties of the 
pressure problem and discuss the necessity of their fulfillment at the discrete 
level. The resulting elliptic problem for the pressure equation is characterized 
by a non-selfadjoint operator. Possibilities of approximate solving the elliptic 
problem are considered using the iterative methods. Special attention is given 
to the numerical algorithms for calculating the pressure on parallel computers. 
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1 Introduction 



> 

l/-v . Mathematical modeling of multicomponent flows in porous media is of great 

I • I importance in oil and gas recovery. Traditionally, the hydrodynamic simula- 

^— s i tors for these applications are based on three-phase black oil model [TlllOj. A 

mathematical model of fluid dynamics in porous media includes differential 
equations, which express the conservation laws of mass and momentum [3l[5]- 
First of all, there are used the continuity equations describing the mass con- 
servation law for each separate phase. The momentum equations in a porous 
r> ' medium are written in the form of Darcy's law, which links the velocity with 
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the pressure. When capillary effects are omitted, the pressure is common to 
all phases. 

Applied mathematical models of mass transfer processes in porous media 
are essentially nonlinear and difficult to study [T61IT9] . Next, it is necessary 
in these models to implement the closure of the system of equations via the 
constant sum of all saturations. Such algebraic components of the model should 
be taken into account in constructing computational algorithms to predict 
multiphase flows in porous media. [HE]. 

Two classes of methods are used to solve approximately unsteady boundary 
value problems for coupled systems of partial differential equations. The first 
of them employs various implicit schemes for the initial system of equations. 
In this case, we face some computational problems in the transition to a new 
time level. The second class of methods reduces computational costs by means 
of using splitting schemes and solving simpler problems at the new time level 
— splitting with respect to physical processes [8lfT3]. The above two classes 
of methods are presented through the fully implicit method (FIM) and the 
implicit pressure explicit saturation (IMPES) approach piHlITU]. 

FIM is widely used in hydrodynamic modeling of oil and gas recovery 
PTj . The fully implicit approximation is used in FIM for all equations of the 
mathematical model. It allows to expect stability of the method and possibility 
to use large time steps. The basic drawback of the method is connected with 
its complexity - we have to solve a large system of nonlinear equations. 

IMPES provides more efficient algorithms for solving the problem at each 
time level. In this approach, we formulate the problem for the pressure with 
implicit approximations in time. After evaluation of the pressure all other un- 
knowns are calculated via explicit approximations. Unfortunately, the problem 
of stability (time step restriction) is typical for the IMPES method. Therefore, 
various modifications of IMPES have been developed in order to improve its 
stability. For example, after evaluation of the pressure we can use implicit ap- 
proximations for the calculation of saturations |20| (the sequential method) . 

In this paper, we highlight the main features of the pressure problem which 
should be taken into account in constructing computational algorithms. There 
are discussed here possibilities of obtaining the pressure equation — elliptic 
for incompressible media and parabolic for compressible ones. If the condition 
of the constant sum of saturations is treated explicitly then the corresponding 
elliptic operator of the pressure problem is non-selfadjoint. This fact should 
be taken into account in constructing iterative algorithms. 

The paper is organized as follows. A basic system of equations is formulated 
in Section [2] to describe multicomponent fluid flows in porous media. This 
mathematical model is obtained at assumptions that capillary and gravity 
forces are negligible. The pressure equation is derived in this section. The 
main features of the grid problem for the pressure are discussed in Section [31 
The simplest uniform grids for the problem in a rectangle are used. 

The emphasis is on iterative methods for calculating the pressure at the 
new time level. The two-dimensional test problem is described in Section [Jj 
The possibility of using standard iterative methods with preconditioners is dis- 
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cussed in Section [S] In sectional we present the results of the iterative solving 
model pressure problems on parallel computers. Conclusions are summarized 
in Section [T] 



2 The pressure problem 

In this section we formulate the basic mathematical model for fluid flows in 
porous media. The system of governing equations for multicomponent flows 
includes the continuity equation for each phase, where a = 1,2, ... ,m — the 
phase index. The mass conservation law for each particular phase is expressed 
by the following equation 

d{<j)baSa) , ■,■ fj^ •, , 10 rT\ 
\-dlV(baUa) ^ -baQa, a = 1,2, ..., TTl. (1) 

Here 4> stands for the porosity, ba is the phase density, Sa — the phase satu- 
ration, Uq — the velocity, and qa — the volumetric mass source. 

For simplicity, we neglect the capillary and gravitational forces. In this 
simplest case the equation of fluid motion in porous media has the form of 
Darcy's law, where the velocity is directly determined by the common pressure: 

k 
Uq = ^ k • gradp, a = 1,2, . . . ,m. (2) 

Ma 

In ^, k is the absolute permeability (in general, symmetric second-rank ten- 
sor) ka — the relative permeability, /!„ — the phase viscosity and p — the 
pressure. 

The unknowns in the system of equations ((ij , ([2]) are the phase saturations 
Sa, a = 1,2, . . . ,m and the pressure (m + 1 unknowns in all). In the simplest 
case, the coefhcients in equations ([1]), ^ are defined as some relations 

(j) = (j){p), ba = ba{p), qa = qa{Sa), K = ka{Sa), IJ-a = COUSt . 

For the sum of saturations of all phases we have 

m 

Y^Sa^l. (3) 

a=l 

After substituting ([2]) in ([T]) and taking into account ([3]), we have a system of 
771 4- 1 equations for to H- 1 unknowns. 

The system of equations H])-© is the basis for the description of multi- 
component flows in porous media. In this system we have not any separate 
equation for the pressure. Equations ([T]) can be considered as the transport 
equation for each phase, whereas the algebraic relation ([3]) can be treated as 
the equation for the pressure. 

Let us consider more convenient forms of system ([Ij-Q, which lead to 
the typical problems of mathematical physics for the pressure. It should be 
noted that such equivalent formulations do exist only at the differential level. 
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At the discrete level such equivalence of formulations is broken even for linear 
problems. So, the choice of the initial form of the equations is essential for 
calculations. 

The most natural way to derive the equation for the pressure is the fol- 
lowing. Divide each equation ^ hy (j>ba > and add them together, which 
gives 




With the natural assumption for compressible fluids 

d{(l)ba) 

> 0, a — l,2,...,m. 

dp 

equation ([4]) for the pressure is the standard parabolic equation of second 
order. In particular, the maximum principle holds for its solutions [7' . 

When using equation Q , the basic system of equations for flows in porous 
media can include m equations 

div k-gradp =-6Qga, (5) 



dt \ fi, 

for Sa, a = 1, 2, . . . , ?n and equation Q for p. In this case equation ([3]) is a 
consequence of (|3]), ([5]). The second approach of common use is connected with 
employing relation ([2]) instead of one of equations ([5]) . For example, equation 
(|4|) is treated as the pressure equation, equations ([5]) are used for Sa, a — 
1, 2, . . . , m — 1 whereas from ^ we get Sm 

m — 1 

Q = l 

Note that the above forms of equations for multicomponent flows in porous 
media are algebraically equivalent only at the differential level. We will try to 
preserve the main points of this equivalence at the discrete level [TH] . 

In the case of variable coefficients cj) ba, the elliptic operator for the pressure 
equation ^ is non-selfadjoint. This fact leads to some problems in using im- 
plicit schemes for equation (j?]). That is why some modifications are employed 
for the pressure equation. For instance, we can obtain the pressure equation 
via the direct summation of equations ([1]) taking into account equation ^ 

E ^^^1^ . f^ div (^^ k . gradp) - f: baqa. (7) 

In this case we have the selfadjoint elliptic operator for the pressure. However, 
this approach has some drawbacks. In particular, this system of equations ([T]), 
(I2|),(IZ1) is not closed, because the basic algebraic relation ^ is not involved. 
Also, we can not say anything about the parabolic property of the pressure 
operator in equation ([T]). 
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3 The properties of the grid operators 

Let us consider stationary and unsteady model problems, which are linear 
prototypes for the pressure problem in modeling multiphase flows. Consider 
the two-dimensional problem in the rectangle 

/? = { X I x= (a;i,X2), 0<xp <lfi, (3 = 1,2}. 

In accordance with Q we solve in f] the boundary problem for the equation 

du 



dt 



^ac(x)£„w = /(x,i), xel2, 0<t<r, (8) 



a=l 



where aQ(x) > £»„, Qa > Q, a — l,2,...,m, and elliptic operators Ca are 
defined by 

under the standard assumptions < Kq, < fc^ < Kq. This equation is supple- 
mented with homogeneous Dirichlet boundary conditions 

w(x,t)=0, xedn, t>0. (10) 

In addition, the initial condition is given in the following form 

u(x,0) =u°(x), xel2. (11) 

In some cases (incompressible media) it is reasonable to consider the sta- 
tionary problem. The boundary value problem is formulated for the equation 

rn 

Au = f{x), ^=^a„(x)/:„, xGf2, (12) 

Q = l 

which is supplemented by the boundary conditions ([TU|) . 

The approximate solution is given at the nodes of the uniform rectangular 
grid in f2: 

(I; = {x I x = (a;i,a;2), xp = iphp, i^ = 0,1,..., iV^, Nphp = lj3, j3 = 1,2} 

and let uj be the set of internal nodes (<D = w U 9a;). For grid functions j/(x) = 
0, X G doj we define the Hilbert space H — L2{uj) with the inner product and 
norm 

{y, w) = Y, y(x)w(x)/ii/i2, ||y|| = (y, y)^/^ 

Approximation in space for problem (I5))- (I11I) will be performed at the 
assumption that the coefficients and solution are sufficient smooth. For the 
elliptic operator Ca we put into the correspondence the grid operator Aa'. 

^"2^ ^ -J^^aixi + 0.5hi, X2)iy{xi + /ii, X2) - j/(a;i, X2)) 
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1 



+T2^"(^i ~ 0-5hi,x2){y{xi,x2) - y{xi - hi,x2)) 



--r2ka{xi,X2 + 0.5h2){y{xi, X2 + /12) - yixi,X2)) 



+ -r^ka{xi,X2 -0.5h2){y{xi,X2) -y{xi,X2-h2)), xGw (13) 
"2 

for all a — 1,2, ...,771. In H the operators Aa, a — 1,2,..., m are selfadjoint 
and positive definite [T^[T5] : 

Aa = A*^, KaSE < Aa < AEaE, a = 1, 2, ..., m, (14) 

where E — the identity operator, and 

d = di + d2, d^^-jsm — — , 
n-/3 2(,3 



4 2 ""^/S 

1-fS 2lp 

The grid operator for the pressure problem can be represented as 



A^Ai + A2, Ap^-jcos'-f, /3=1,2. 
hi 2lr-. 



A^'^aa{Ti)Aa, X e w. (15) 

a=l 

In general (non-constant coefficients aQ.(x), a — 1,2, ...,m) the operator A is 
non-selfadjoint. It approximates the corresponding differential operator with 
the error of 0(|/ip), where |ft.p ^ hj + hj. 

After discretization in space we go from (|5|)-([TT|) to the differential-operator 
equation 

^+Ay = f{t), 0<t<T, (16) 

considered on the set of grid functions y{t) G H. The initial condition is taken 
in the form 

y(0) = u°. (17) 

For the stationary problem ([TU]). (IT^ the grid analog has the form 

Ay = /. (18) 

To solve approximately problem (IT5|) . (IT7|) . we use the standard two- level 
schemes. Let r be the fixed time step and y" = y(i"), i" = ?^t, n = 0, 1, ..., A^, Nt 
T. Equation (J16D is approximated by the two-level scheme with weights 

n+l _ n 

^ ^ + A(a2;"+i + (l-a)y") = ^", n = 0, 1, ..., A^ - 1, (19) 

T 

where, for example, ip" = fiaf"^^ -|- (1 — (T)t"). It is supplemented by the 
initial condition 

y" = u". (20) 
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The difference scheme (|T^ . (EU]) has the approximation error in time 0{t^ + 
(cr-0.5)r). 

If we employ the fully implicit scheme (u = 1), then the transition to the 
new time level is performed through solving the grid problem 

^E + A^y^f. (21) 

The main subject of our consideration is the methods of solving grid problems 
(|18p and (|2ip . which are linear prototypes for stationary and unsteady problem 
for the pressure. The primary question here is the non-selfadjoint property of 
the grid operator A. 

For the grid problem (fTSj) the maximum principle holds [T2] . With regard to 
considered approximations on the five-point stencil, we formulate it as follows 
[l4] . Consider the difference equation 

7(x)y(x) -ai(x)y(a;i - hi,X2) - I3i{^)y{xi +/ii,X2)- 

- a2(x)y(a;i,a;2 - ft.2) -/32(x)2/(xi, 2:2 -1-/12) = </5(x), x G w, (22) 

which is supplemented by boundary conditions 

y(x) = 0, X e duj. (23) 

We assume, that the coefficients of the difference scheme (P^ satisfy the con- 
ditions 

aj(x) > 0, /3j(x) > 0, J = 1, 2, 7(x) > 0, x G c^. (24) 

Let in the difference scheme (P^ - (P5|) we have f{x.) > for all x G w (or 
(^(x) < for X G w). Then for 

7(x) >ai(x)-f a2(x)-H/3i(x)+/32(x), xGw (25) 

we have (the grid maximum principle) 2/(x) > for all x G a; (2/(x) < for 
X G w). In our case (see p^ . (fTS]) ) fulfillment of the sufficient conditions (|^ 
can be verified directly. Because of this, for the grid operator at the new time 
level (PT|) we have the strict diagonal dominance. 

To study properties of operators A and A in Hilbert spaces Ti. = L2{fi) and 
H = ^2(0;), it is convenient to treat A and A as the corresponding convection- 
diffusion operators. In this case it is possible to employ in our research the 
results from [9l[T4]. 

Taking into consideration ([9|) and (fT2)) . we have the representation 

m 

A=Y,-^<^' Aa=T>a+Ca, a = l,2,...,m, (26) 



where 



VaU = —div{da{x)gTa.du), (27) 

CaU — Wq gradu. (28) 
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The effective diffusion coefScient and convection velocity for the separate phase 
a are 

Then the pressure operator takes the form of convection-diffusion operator 
with the convective term in the non-divergent form. Note that apphcation of 
equation ([7]) to evaluate the pressure corresponds to using only the diffusion 
part (gZl) of the operator ^. 

Operators of diffusion in the above assumptions about the coefficients are 
self-adjoint and positive definite in T-L = L2{fi). Next, we present some facts 
about the properties of convective transport operators. A detailed discussion 
of these issues is given in the book [14] . 

We have the following representation 

Ca = Cq - -divwa, a^l,2,...,m, (29) 

where Cq is the operator of convective transport in the symmetric form: 

CaU — - (wq gradu + div(wQ,u)) . (30) 

The operator Ca is skew-symmetric in H: 

C„ = -C* (31) 

for any "Wq,(x), x £ i7. 

From (l^nj) and (I5T|) we directly obtain the estimate for the energy of the 
convective transport operator Ca'- 

\{C^u,u)\<Ma\\u\\^ (32) 

Ma^ -^\\divwa\\c{n), l|w||c(fi) = max |ii(x)|. (33) 

It is interesting to consider the subordination estimate for the operator of 
convective transport with respect to the diffusion operator. In our model two- 
dimensional problem the corresponding estimate has the following form 

\\CaU\\^ <Ma{Vc,U,u), (34) 

at Wc = (wa , Wa ) with constant 



2 

Ada < max ■ 

QaKa /3=1,2 



{wi^')' 



(35) 
an) ) 



These properties of differential operators of diffusion and convection (j29p , 
dnU, (1321), (EH) are inherited not only for the difference operators on rectan- 
gular grids [14], but also for difference operators on irregular grids with the 
Delaunay triangulation |17j . We consider this issue here for the grid operator 
(fnjl . (dU). First of all, we are interested in the grid analog of (l^ - (P5|) . 
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Taking into account 

-a(^)i (k{x + h/2)y±±4^y^ - Hx h^y^''^ ~ ^^^ - '^^ 



h \ h h 

1 f a{x + h) + a{x) y{x + h) - y{x) 
—h \ 2 ^(" + ^'^^ h 

a{x) + a{x ~ h) ,^^ ui^\ y^^) -yjx-h) 
_ k{x - n/l) - 

- - k{x + h/2) -^ + 

1 a{x) - a{x - h) ,^^ y{x)--y(^c-h) 

2 h ''^'^ ~ ^'^^ h ' 
similarly (1261) we obtain 

A=^Aa, Aa=Da_^Ca, a = l,2,...,m. (36) 

The grid diffusion operator has the form 

aa{xi+h-Y,X2) + aa{xi,X2) . , ^ c f, M I , U \ i \\ 

ka{xi + 0.bhi,X2)[y(xi + hi,X2) -y{xi,X2)) 

ka{xi - 0.5hi,X2){yixi,X2) -y{xi - hi,X2)) 

ka{xi,x2 + 0.5h2){y{xi,x2 + /12) -y{xi,x2)) 



H '■ TTa ' kaixi,X2 - 0.5h2){y{xi,X2) - yixi,X2 - /l2))- 

(37) 
Similarly (fH)) in H = L2{lo) we have 

Do. = Dl, D^> p^K^SE. (38) 

Approximation of the convective part of the grid operator A is conducted 
via setting the coefhcients (effective velocity w^) on the grids shifted in the 
corresponding direction on the half-step. Let us define with accuracy of Od/ip) 
the components of the grid analog of Wq using the following relations 

(1). , „ ^, . aaixi+hi,X2)-aa{xi,X2) , . 

w^^>[xi +0.5hi,X2) = r ka(xi +0.5hi,X2), 

n-i 

(2)/ , nKh\ aa{xi,X2 +h2) -aa(xi,X2) ,^ , , nru \ (nn\ 

w)^'[xi,X2 +O.5/12) = r ka{Xi,X2 + 0.5/l2)- (39) 





2hl 






aa{xi, 


.X2) + aa{xi 


~h,. 


,X2) 




2hl 






aa{xi, 


,X2 + h2) + a 


q(xi 


,X2) 




2hl 






aa{xi 


,X2) +aa{xi, 


X2 - 


112) 
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The convective transport operator in representation p6p has the form 
CaV = TT^a (2^1 + 0.5/11, X2) r \- 

Z 111 

1 (1)/ r,Kh ^ yi^u^2) -yjxi -hi,X2) 
-w^^(a;i -0.5/11,2:2) r ^ 

Z /ll 

-W^ ^(Xi, 2:2 + 0.5/12) r h 

Z /l2 

1 (2)/ nr;^ ^ ^(a^l, X2) - ^(xi , 2:2 - ^2) ..„^ 

-u;^^(xi, 2:2 -0.5/12) r • (40) 

Z 112 

The grid analogue of ^TE\\ can be written as 



Ca^Ca- -div,iW„, 


(41) 


where 

Wq (xi +0. 5/11,^2) - Wq (2:1 - 0.5/ii,X2) 




dlA- h Wq, — 1 




Wa (a;i + ,X20.5/i2) - Wq^(xi,2;2 - O.5/12) 


(42) 


/l2 


For the skew-symmetric part 




L'Q = ~^a 


(43) 



we have 

CaV = 777— wi^^a^i +0.5hi,X2)y{xi +hi,X2)- 
Zfli 

— — w£^^(xi - 0.5/ii,X2)y(xi - /ii,X2)-|- 
z/li 

——W^^^ (2:1, 2:2 + 0.5/12)2/(^1, X2 +/i2)- 
Z/I2 

— -u;(f)(xi, X2 - 0.5h2)yixi,X2 ~ /i2)- (44) 

The following grid analog of (l32t takes place: 

|(ay,2/)|<M„||yf, (45) 



where now (see (1331) ') 

^^a = oll^^^'''*^"llc(i.^)> Il2/I!c(cu) =max|y(x)|. (46) 

The subordination inequality (see ([Mj) and ([55l) ) has the form 

||(C„yf <M„(i?a2/,y), (47) 
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C(u]) C(u 

(48) 
The fundamental issue here is that for these approximations the constants Ma 
and Ma are the complete grid analogues of the corresponding constants Ada 
and Aia for the differential problem. 



4 The test problem 

Capabilities of iterative methods for approximate solving the pressure equation 
in modeling multiphase flows in porous media are illustrated here using the test 
grid problem. We consider equation (PT|) . which corresponds to the calculation 
of one time step in the numerical solution of problem ([5])- pT|) . Numerical 
experiments are conducted for problem (j2ip with / = 1 in the unit square 
[Ifi = 1, /3 = 1, 2) on the grid h = hi = h2 {N = Ni = N2). 

Particular attention should be given to the coefficients of equations ([5]) , (jH]) 
in order to take into account peculiarities of these problems, namely, inhomo- 
geneity of aQ(x), a = 1, 2, ..., m. Taking into account Q, we set 

fca(x) - — -— , a == 1,2, ...,m. 
a«(x) 

Consider two-phase medium (m = 2) with an incompressible fluid as the first 
phase 

ai(x) = 1, fci(x) = 1. 

Compressibility of the second phase is defined as follows: 

a2(x) = exp(-e((.T^i - 0.5)2 ^ ^^^ _ q 5^2)^^ 
fc2(x) - r/exp(e((xi - 0.5)2 ^ (^^ _ o_5)2))_ 

The diffusion part of operator (P7|) is 

Piw = — divgradw, 232^ = ^'7divgradu. 

Properties of the considered problems are defined(see (UHl)) by the vectors 
Wq, a = 1, 2, ..., m. For the test problem we have 

wi =0, W2 = {-2ri^{xi - 0.5), -2ri^{x2 - 0.5)). 

In this case 

div W2 = — 477^. 
For the constants in the estimates ([5^ and p4p we obtain 

M2 = 2?7|ei, M2 = 2r7C'- 

Thus, the governing numerical parameters for this problem are 77 and £,. The 
sign of ^ can be any, moreover, it defines the fundamental difference in the 
behavior of the solution (the pressure) in the vicinity of the production or 
injection well. 
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5 Iterative solution of the problem 

For numerical solving the test problem we use iterative methods. In the cor- 
responding grid equation (j21l) the operator A is non-selfadjoint. Therefore, we 
use iterative methods for grid problems with unsymmetric matrices 111115). 
The standard Generalized Minimal Residual Method (GMRES) with different 
preconditioners has been employed. 

To solve the test problem, the PETSc library [5] has been used. The PETSc 
library, developed in the Argonne National Laboratory, is a powerful set of 
freely available multi-platform compatible tools for the solution of large-scale 
problems governed by partial differential equations. Experiments were carried 
out with the following preconditioners: 

none — without preconditioning; 

jacobi — the Jacobi method; 

sor — the successive overrelaxation method; 

ilu — the incomplete LU factorization; 

mg — the multigrid method. 

Table [1] shows the dependence of the computational cost (the number of 
iterations) on the physical parameters of the problem. Features of the problem 
are clearly defined by the parameters 77 and ^. The calculations were performed 
using the unpreconditioned GMRES method on the grid with 256 x 256 un- 
knowns. We see that with increasing of r/ and/or |^| the number of iterations 
decreases. The same is true for negative values of rj. 



Table 1 The dependence of the number of iterations on the physical parameters (g, rj) 



V 


-10 


-1 


0.0 


1 


10 


0.01 


4710 


4686 


4682 


4678 


4629 


0.1 


4413 


4726 


4699 


4652 


3828 


1 


1445 


4790 


4790 


4289 


1689 


10 


879 


4568 


4884 


3919 


1088 


100 


857 


4493 


4903 


3856 


1026 



Effect of preconditioning on different grids is shown in Table[2j Calculations 
were performed at 77 = 1. It is easy to see that the multigrid preconditioner is 
the best. 

In addition, it is interesting to look at the effect of the time step r. The 
unpreconditioned GMRES method was used with the grid of 256 x 256 un- 
knowns. From Table [3] we see that the number of iterations decreases with r. 
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Table 2 The number of iterations for different 


preconditioners depending 


on the grid 


grid 


preconditioncr 


-10 


-1 


0.0 


1 


10 




none 


538 


1283 


1270 


1150 


523 




jaeobi 


507 


1283 


1270 


1150 


517 


128 X 128 


sor 


222 


280 


284 


281 


156 




ilu 


175 


217 


214 


215 


128 




mg 


5 


5 


5 


5 


5 




none 


1445 


4790 


4790 


4289 


1689 




jaeobi 


1443 


4789 


4790 


4284 


1675 


256 X 256 


sor 


350 


807 


765 


703 


389 




ilu 


325 


609 


566 


534 


294 




mg 


5 


5 


5 


5 


5 




none 


3271 


17685 


18777 


16890 


6172 




jaeobi 


3429 


17721 


18777 


16873 


6105 


512 X 512 


sor 


1043 


2699 


2987 


2510 


1120 




ilu 


764 


2050 


2045 


1596 


828 




mg 


5 


5 


5 


5 


5 



Table 3 The immber of iterations for various r 



1] 


T 


-10 


-1 


0.0 


1 


10 




0.01 


864 


866 


866 


866 


868 




0.1 


3329 


3324 


3323 


3321 


3306 


0.01 


1 


4710 


4686 


4682 


4678 


4629 




10 


4913 


4887 


4882 


4877 


4822 




100 


4936 


4907 


4903 


4899 


4844 




0.01 


892 


924 


926 


928 


933 




0.1 


3203 


3420 


3412 


3396 


3015 


0.1 


1 


4413 


4726 


4689 


4652 


3828 




10 


4568 


4914 


4884 


4830 


3919 




100 


4687 


4934 


4903 


4850 


3939 




0.01 


958 


1413 


1444 


1446 


1085 




0.1 


1313 


3927 


3953 


3661 


1614 


1 


1 


1445 


4790 


4790 


4289 


1689 




10 


1455 


4896 


4893 


4368 


1695 




100 


1456 


4908 


4904 


4374 


1696 



6 Parallel implementation 

The parallel formulation is based on the domain decomposition methods. The 
main idea of these methods is to divide the original computational domain into 
subdomains. A separate processor, which is identified by its rank, is assigned 
to each subdomain in order to perform the computations. For inter processor 
communications the Message Passing Interface (MPI) is used. 

The systems of linear equations are solved by the parallel version of the 
preconditioned GMRES algorithm. In our computations the none, bjacobi 
(doing the ILU-factorization of a local part of the matrix at each processor) 
and multigrid preconditioners were used. The calculations were performed on 
the grid 512 x 512 at 77 = 1. 
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The parallel code was run on a cluster of North-Eastern Federal University. 
The cluster consists of four computing nodes, each node has two quad-core 
processors Intel Xeon E5450 (3.00 GHz) with 16 Gb RAM. 

The results of the parallelization efficiency of computations are given in 
Table 21 The table shows the estimation of computational costs, since the 
number of iterations is almost independent of the number of running processes. 

Table 4 Computation time in seconds, np — number of processes, pc — preconditioner 



np 


pc 


-10 


-1 





1 


10 


16 


none 

bjacobi 

mg 


41.26 
11.48 
3.06 


226.68 
31.16 
3.17 


241.78 
31.69 
3.10 


202.67 
24.94 
3.10 


72.21 
12.66 
3.14 


8 


none 

bjacobi 

mg 


57.35 
13.55 
2.38 


304.48 
43.54 
2.39 


236.72 
28.92 
3.64 


294.56 
3502 
2.85 


108.07 
12.46 
2.91 


4 


none 

bjacobi 

mg 


86.28 
29.08 
4.99 


593.00 
86.85 
4.15 


583.87 
88.96 
3.38 


538.89 
69.80 

4.84 


220.66 
34.30 
4.92 


2 


none 

bjacobi 

mg 


216.75 

75.55 

5.9 


1328.398 
177.12 

5.87 


1361.5 

145.78 
7.98 


883.42 

146.34 

5.90 


471.60 
79.90 
7.95 


1 


none 

bjacobi 

mg 


315.07 
73.33 
9.23 


1686.16 

197.30 

8.55 


1799.68 
216.68 

8.57 


1612.7 

200.49 

8.56 


590.24 
86.07 
8.54 



7 Conclusions 

1 . The basic features of the pressure problem associated with the non-selfadjoint 
operator are considered for multiphase flows in porous media. 

2. It was found that the computational cost of solving the model pressure 
problem does not depend strongly on ^, more pronounced dependence is 
on the physical parameter 77. This means that the number of iterations 
depends basically on the various properties of the phases than on the value 
of external sources. 

3. Parallel computations have been performed using standard techniques with 
various preconditioners. 
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